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Abstract 

The results obtained from molecular dynamics simulations of the friction at an interface between 
polymer melts and weakly attractive crystalline surfaces are reported. We consider a coarse-grained 
bead-spring model of linear chains with adjustable intrinsic stiffness. The structure and relaxation 
dynamics of polymer chains near interfaces are quantified by the radius of gyration and decay of 
the time autocorrelation function of the first normal mode. We found that the friction coefficient 
at small slip velocities exhibits a distinct maximum which appears due to shear-induced alignment 
of semifiexible chain segments in contact with solid walls. At large slip velocities the decay of the 
friction coefficient is independent of the chain stiffness. The data for the friction coefficient and 
shear viscosity are used to elucidate main trends in the nonlinear shear rate dependence of the slip 
length. The influence of chain stiffness on the relationship between the friction coefficient and the 
structure factor in the first fluid layer is discussed. 

PACS numbers: 68.08.-p, 83.80.Sg, 83.50.Rp, 47.61.-k, 83.10.Rs 
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I. INTRODUCTION 



Understanding the interfacial rheology of complex fluids is important in many processes 
relevant to technological applications including polymer processing 2!, boundary lubri- 
cation ^, and dewetting of polymer films |4|. Numerous experimental studies have demon- 
strated that flow velocity profiles in nanoconfined systems can be significantly influenced by 
slip at polymer-solid interfaces js, 6|. The measure of slippage is the so-called slip length, 
which is defined as the distance between the physical interface and imaginary plane where 
the extrapolated velocity profile reaches the substrate velocity. During the last two decades, 
the dependence of the slip length on flow conditions and material properties of substrates 
was extensively investigated by molecular dynamics (MD) simulations for monatomic 



and polymeric 



7H26| 



27|-|45| fluids. 



In the case of simple shear flow illustrated schematically in Fig.[Tl the shear stress is the 
same in the bulk of the channel and at the interface; and, therefore, the slip length can 
be calculated from the ratio of the fluid viscosity to the friction coefficient at the interface 
according to 

L. = f , (1) 

where the friction coefficient k is defined by the relation between the slip velocity and wall 
shear stress At equilibrium, both the fluid viscosity 



48| and friction coefficient [1^, 



16 



24 



49| can be estimated from the Green-Kubo relations, and the slip length in the limit 



of zero shear rate is then computed from Eq. ([T]). In the presence of flow, the fluid viscosity 
and friction coefficient might depend on shear rate and slip velocity respectively: and, as a 



result, the slip length is often a nonlinear function of shear rate 10 



14 



resp 
22-I24, 



36, 



38, 



42i. For 



example, it was recently shown that at an interface between unentangled polymer melts and 
passive surfaces (no chemical bonds with the surface and weak wall-fluid interaction energy), 
the slip length passes through a local minimum at low shear rates and then increases rapidly 



at higher shear rates 36|, |42]. This non-monotonic behavior was explained by computing 



the rate-dependent viscosity and the friction coefficient that undergoes a transition from a 



36, 



42|. One of the 



constant value to the power-law decay as a function of the slip velocity 
motivations of the present study is to examine whether these conclusions hold for different 
simulation ensembles and intramolecular potentials. 

In the past two decades, a number of MD studies have demonstrated that the degree 



2 



of slip at the interface between molecular liquids and crystalline surfaces depends on the 



structure of the first fluid layer in contact with the periodic surface potential 



36 



39 



0, Qjie 



42| . In general, the interfacial slip is suppressed with increasing height of the peak 



in the fluid structure factor computed at the main reciprocal lattice vector. In turn, the 
in-plane order within the adjacent fluid layer is determined by several factors, including 



the commensurability o: 



energy 



B Q, Q, Q, 



42 



43j, and fluid pressure 



11 



28 



39 



42|, wall-fluic 



the liquid and solid structures 



16 



interaction 



28 



35 



42 



42| . Interestingly, recent simulation results have shown 



that the friction coefficient in the linear slip regime is a function of a combined variable 
that is a product of the height of the main peak in the structure factor and the contact 
density of the first fluid layer near the solid wall j^]. However, at present, there exists no 
exact relationship between the friction coefficient (or the slip length) and the microscopic 
properties of the liquid-solid interface. 

The formation of the interfacial fluid layer and hydrodynamic boundary conditions can 
be significantly affected by the polymer chain architecture. For instance, several MD studies 
have reported flow profiles with a finite slip velocity in thin alkane films where stiff polymer 



chains tend to align in layers parallel to the surface 



31 



3311. On the other hand, a weaker 



density layering near the wall and pronounced slippage were observed for branched alkane 



molecules 



34| . Notably, it was demonstrated that upon increasing the length of linear, freely- 



jointed chains, the structure of the first fluid layer near a solid wall is reduced, resulting in 
smaller values of the friction coefficient and larger slip lengths at low shear rates [35]. Later, 
it was found that the existence and location of a double bond along the backbone of linear 
oligomers affect the structure of the interfacial fluid layer near crystalline aluminum walls 
and lead to either negative or large positive slip lengths More recently, it was shown 
that slip velocity is reduced for liquids which consist of molecules that can easily conform 
their atoms into low-energy sites of the substrate potential 43|]. Despite extensive research 
on the fluid structure and shear response in thin polymer fllms, it is often difficult to predict 
even qualitatively the influence of liquid molecular structure on the interfacial slip. 

In this paper, we investigate the effect of chain bending stiffness on the fluid structure 
and friction coefficient at interfaces between linear polymers and crystalline surfaces. We 
will show that the relaxation dynamics near weakly attractive surfaces is significantly slowed 
down for stiffer chains at equilibrium. The orientation of semiflexible chains in shear flow 
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leads to the enhanced density layering away from the walls and partial alignment of extended 
chain segments in the first fluid layer. It will be demonstrated that the shear-induced 
ordering of the chain segments produces a distinct maximum in the liquid structure factor 
and the friction coefficient at small slip velocities. Finally, the characteristic features in 
the shear rate dependence of the slip length are interpreted in terms of the shear-thinning 
viscosity and the dynamic friction coefficient. 

The rest of this paper is organized as follows. The details of molecular dynamics simu- 
lations, interaction potentials, and equilibration procedure are described in Section |Tll The 
simulation results are presented in Section IIIII More specifically, the chain conformation and 
relaxation dynamics are analyzed in llll A[ examples of density and velocity profiles are pre- 
sented in llll B[ shear viscosity and slip lengths are reported in llll C[ the velocity dependence 
of the friction coefficient is examined in IIIIDl and, finally, the fluid structure near solid walls 
is considered in llll E[ Brief conclusions are given in the last section. 



II. MOLECULAR DYNAMICS SIMULATION MODEL 



We consider a coarse-grained bead-spring model of unentangled polymer melt, which 
consists of M = 480 linear chains of = 20 beads (or monomers) each. In this model any 
two fluid monomers interact via the truncated Lennard- Jones (LJ) potential 



VlAt) 



Ae 



r . 



r . 



(2) 



where e and a are the energy and length scales of the fluid phase. The cutoff radius Vc = 2.5 a 
and the total number of fluid monomers Nf = 9600 are fixed throughout all simulations. 
Similarly, fluid monomers interact with wall atoms via the LJ potential with the following 
parameters e^f = 0.8 e, o"„f = a, and = 2.5 cr. 

Any two consecutive beads in a polymer chain interact through the finitely extensible 
nonlinear elastic (FENE) potential 50 1 



VpENEir) 



(3) 



5l| . The combination of LJ and 



with the standard parameters = 30ea~^ and Vo = 1.5 a 
FENE potentials yields an effective bond potential between the nearest-neighbor beads with 
the average bond length b = 0.97 a 5l|. This bond potential is strong enough to prevent 
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chain crossing and breaking even at the highest shear rates considered in the present study. 
In addition, the flexibihty of polymer chains is controlled by the bending potential as follows: 



UbendiO) = kg {1- cos 9), (4) 

where kg is the bending stiffness coefficient and 6 is the angle between two consecutive bonds 



10 



28 



54 



55|. 



along a polymer chain [52|]. In the present study, the bending stiffness coefficient was varied 
in the range ^ kg ^ 3.5 e. A snapshot of the confined polymer melt that consists of 
semifiexible linear chains with the bending coefficient kg = 2.5 e is shown in Figured 

In order to remove viscous heating generated in the shear flow, the motion of fluid 
monomers was coupled to an external heat bath via a Langevin thermostat applied in 
the y direction to avoid bias in the shear flow direction (the x direction). This is a standard 
thermostatting procedure often used in MD simulations of sheared fluids 
Thus, the equations of motion for fluid monomers are summarized as follows: 

-^^ = -E^' (5) 

myi + mViji = - — ^ + fi , (6) 
^ dyi 

where F = 1.0 r^^ is the friction coefficient that controls the damping term, Vij is the total 
interaction potential, and /j is a random force with zero mean and variance {fi{0)fj{t)) = 
2mkBTT5{t)5ij obtained from the fluctuation-dissipation theorem. The Langevin thermo- 
stat temperature is set T = Ller/fc^, where ks refers to the Boltzmann constant. The 
equations of motion were solved numerically using the fifth-order Gear predictor-corrector 



algorithm 



56 1 with a time step At = 0.005 r, where r = ^Jma"^ /e is the LJ time. Typ- 



ical values of the length, energy, and time scales for hydrocarbon chains are a = 0.5 nm. 



£ = 30 meV, and r = 3 x 10"^^ s 



56|. 



The polymer melt is confined between two crystalline walls as illustrated in FigureO 
Each wall consists of 1152 atoms distributed between two layers of the face-centered cubic 
(fee) lattice with density = 1.40 cx"^. For computational efficiency, the wall atoms are 
fixed rigidly to the wall lattice sites, which form two (111) planes with [112] orientation 
parallel to the x direction. The nearest-neighbor distance between the lattice sites within 



5 



the (111) plane is d = 1.0 a and the first reciprocal lattice vector in the x direction is Gi = 
(7.23 a~^, 0). The channel dimensions in the xy plane are measured to be = 20.86 a and 
Ly = 24.08 cr. Periodic boundary conditions were applied along the the x and y directions 
parallel to the solid walls. 

In our simulations, the distance between the wall lattice planes, which are in contact 
with the fluid phase, was flxed at h = 22.02 a. Hence, the volume accessible to the fluid 
phase corresponds to the fluid monomer density p = Nf/LxLy{h — a) = 0.91 a"^; and, in the 
absence of shear flow, the resulting fluid pressure and temperature are 1.0 e cr^^ and 1.1 e/kB 
respectively. In the present study, the relatively low polymer density (or normal pressure) 
was chosen based on the results from our previous study where it was shown that for weak 
wall-fluid interactions and p ^ 1.02 (or P ^ 5.0 ea^^), the fluid velocity proflles remain 



linear in a wide range of shear rates |36|. In contrast, it was demonstrated that at higher 
polymer densities (or pressures), the velocity profiles acquire a pronounced curvature near 
the wall and the relaxation of flexible polymer chains in the interfacial region becomes very 
slow 



The system was flrst equilibrated for about 5 x lO^r while both walls were at rest. Then, 
the velocity of the upper wall was increased gradually up to a target value, followed by an 
additional equilibration period of about 5 x lO^r. In this study, the upper wall velocity was 
varied over about three orders of magnitude 0.005 ^ Ur/a ^ 5.5. Once the steady shear 
flow was generated, the velocity, density, and temperature proflles were averaged within 
horizontal bins of thickness Az = 0.01 a for a time period up to 5 x lO^r. At the lowest 
upper wall speed, U = 0.005 a/r, the velocity proflles were computed in 24 independent 
systems for the total time period of about 5 x lO^r. An upper estimate of the Reynolds 
number at high shear rates is Re = phU/fi = 0(10), which is indicative of laminar flow 
conditions in the channel. 

We flnally note that MD simulations were also performed at a constant normal load, 
where the distance between the walls was allowed to vary under the constant normal pressure 
Pj_ = l.Oe cr~^ applied to the upper wall. However, we did not observe any qualitatively new 
behavior; and for the sake of brevity, these results are not reported in the present study. 
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III. RESULTS 



A. Chain conformation and relaxation dynamics 

The spatial configuration of polymer chains is well characterized by the radius of gyration, 
which is defined as the average distance between monomers in a polymer chain and its center 
of mass as follows: 



1 ^ 

i=l 



where rj is the position vector of the i-th monomer, = 20 is the number of monomers per 
chain, and rem is the chain center of mass defined as 

1 ^ 

i=l 

Figure [3] shows the radius of gyration and its components along the ct, y, and z directions 
as a function of the bending stiffness coefficient. The chain statistics were collected in the 
interfacial regions (where at least one monomer in a chain is in contact with wall atoms) 
and in the middle of the channel (the chain center of mass is located further than 6 a away 
from the walls). As expected, in both cases Rg increases with increasing chain stiffness. 
Note that even at the largest value of the bending stiffness coefficient = 3.5 e, the size 
of polymer chains is smaller than the channel dimensions. The simulation results in Fig. [3] 
indicate that in the bulk region the chain configuration is isotropic, while near the interfaces 
polymer chains become flattened, i.e., Rgz < Rgx ~ Rgy, which is in agreement with previous 



MD studies 57H59|. It is also apparent that fully flexible chains in contact with the walls 
are packed on average within the first two fluid layers [2Rgz ~ 2 a for fcg = in Fig.[3](a)]. 
In contrast, semiflexible chains extend up to about three molecular diameters from the 
walls. Visual inspection of the polymer chains in the interfacial regions revealed that the 
conformation of semiflexible chains consists of locally extended segments within the first fluid 
layer and segments of several monomers oriented away from the walls. Finally, regardless of 
the chain stiffness, the total radius of gyration is nearly the same in the bulk and close to 
the walls due to the relatively weak wall-fluid interaction energy. 

The local relaxation dynamics in confined polymer films can be described by the decay 



of the time autocorrelation function of normal modes 



39 



58 



60 



6l| . By definition, the 
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normal coordinates for a polymer chain that consists of monomers are given by 

x.w = ^Er^w-^^^r^' (10) 

i=l 

where Tj is the position vector of the z-th monomer in the chain, and p = 0,1, N — 1 is 
the mode number 6^. The longest relaxation time corresponds to the first mode p = 1, 



i.e., to the relaxation of the whole chain 



62| . The normalized time autocorrelation function 



for the first normal mode is then defined as follows: 

C,{t) = (Xi(t) ■ Xi(0))/(Xi(0) ■ Xi(0)). (11) 

In our study, the autocorrelation function [Eq. ( fTTI) ] was computed separately in the inter- 
facial regions (where the chain center of mass is confined within 3 a from the walls) and in 
the bulk of the channel where the fluid density is uniform (the center of mass is located at 
least 6(7 away from the walls). An important aspect is that the autocorrelation function 
was averaged only for those polymer chains whose centers of mass remained within either 
the interfacial or bulk regions during the relaxation time interval. 

FigureH] shows the relaxation of the time autocorrelation function at equilibrium (i.e., 
when both walls are at rest) for selected values of the bending stiffness coefficient. As is 
evident from Fig.|l](a), the relaxation rate of polymer chains in the bulk region decreases 
with increasing bending stiffness. A similar effect was reported previously for linear bead- 
spring chains with variable bending rigidity {gsI, indicating that the reorientation dynamics 
in the melt is slowed down for more rigid polymer chains. In our study, the decay rate 
of the autocorrelation function in the bulk is well described by the exponential function 
Ci(t) = exp(— t/ri), where ri is the characteristic relaxation time. The inset in Fig.|l](a) 
presents the variation of ri as a function of the bending stiffness coefficient. The inverse 
relaxation time, I/ti, is related to the characteristic shear rate, above which the shear 
viscosity is expected to exhibit non- Newtonian behavior (see discussion below). We note that 
for stiffer chains, the estimated shear rate is about 2.5 x 10~^r~^, which is about the lowest 
shear rate accessible in coarse-grained MD simulations (without excessive computational 
time requirements). 

In contrast, the decay in time of the autocorrelation function is much slower for semiflex- 
ible polymer chains in the interfacial regions, see Fig.|4](b). Similar results were observed 
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previously in MD simulations of freely-jointed 5-mers adsorbed on weakly physisorbing sur- 
faces, i.e., the relaxation time of adsorbed chains is about an order of magnitude larger than 
in the bulk 60|. In our setup, the relaxation time of flexible chains in the interfacial region 



is only slightly larger than in the bulk, which means that the relaxation dynamics is weakly 
affected by the substrate. With increasing chain stiffness, however, the rotational relaxation 
is significantly slowed down. The data shown in Fig.|l](b) cannot be well fitted by the single 
exponential function. We comment that the typical relaxation time for polymer chains in 
the interfacial regions was used to determine the time interval for averaging the radius of 
gyration and bond orientation. It should also be mentioned that test simulations of polymer 
chains with larger stiffness coefficients, kg = 4.0 e and 4.5 e, have shown that their relaxation 
dynamics near interfaces is extremely slow and cannot be accurately resolved (results not 
reported). 



B. Fluid density and velocity profiles 

The averaged monomer density profiles are presented in Fig.[5]for fully flexible {kg = 0.0 e) 
and semiflexible {kg = 3.0 e) polymer chains at small {U = 0.01 cr/r) and large {U = 4.0 a/r) 
upper wall velocities. Near the solid walls, these proflles exhibit typical density oscillations 
that gradually decay to a uniform proflle in the middle of the channel. Notice that the 
magnitude of the flrst peak in the density proflles (deflned as the contact density) is higher 
for stiffer polymer chains. For example, the contact density is pc = 3.19 cr^'^ for flexible 
chains and pc = 3.44 for kg = 3.0 e when the upper wall velocity is U = 0.01 a/r in 
Fig.O This result, at flrst glance, appears to be somewhat counterintuitive because one 
might expect that flexible chains can pack more effectively near a flat surface. However, 
with increasing bending rigidity, the persistence length of polymer chains increases; and, 
therefore, the flrst fluid layer contains more extended chain segments. When U = 0.01 cr/r 
in Fig.[5l the average number of consecutive monomers per polymer chain in the flrst fluid 
layer is Ng^g ~ 3.1 for flexible chains and Ngeg ~ 5.2 for kg = 3.0 e. It turns out that these 
locally extended chain segments arrange themselves more tightly near the surface. We also 
note that similar trends in the fluid density layering were reported in other coarse-grained 
MD simulations; namely, that with increasing length of (semi)flexible polymer chains, the 
amplitude of density oscillations near a solid wall becomes (larger) smaller {32! . IsT^ . 
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As shown in Fig.[5](a), the height of the density peaks in the case of flexible chains is 
reduced at the higher upper wall speed U = 4.0 a/r. At these flow conditions, the slip 
velocity of the first fluid layer is relatively large (of about 1.0 a/r), and the temperature of 
the fluid near the walls is higher than the temperature of the Langevin thermostat, leading 
to a reduced density layering. This is consistent with the results of previous MD studies 
where the shear response of thin polymer films was examined in a wide range of shear 



rates 



36 



42| . Interestingly, while the amplitude of the first two peaks in the density 



profile for semiflexible chains {ke = 3.0 e) is also reduced at the higher upper wall speed 
U = A.Oa/r, the orientation of more rigid chain segments along the shear flow direction 
produces slightly higher density in the 3rd, 4th, and 5th fluid layers [see Fig.[n](b)]. 

The representative velocity profiles are plotted in Fig. for the lowest U = 0.005 cr/r 
and intermediate U = 0.5 cr/r upper wall speeds and kg = 0.0 e, 2.0 e, and 3.0 e. For 
U = 0.005 cr/r, despite extensive averaging, the data remain noisy because the average flow 
velocity is much smaller than the thermal fluid velocity vt = kBT/m. In all cases, the 
velocity profiles are anti-symmetric with respect to the center of the channel and linear 
except within about 2 a near the walls. Surprisingly, the dependence of slip velocity on 
bending stiffness shows opposite trends for the reported upper wall speeds; namely, the slip 
velocity for flexible chains is smaller for U = 0.005 cr/r in Fig.[6](a), while it is larger for 
U = 0.5 cr/r in Fig.|6](b). This result illustrates that the effect of chain stiffness on the 
interfacial slip strongly depends on flow conditions. In what follows, the shear rate was 
extracted from the linear part of velocity profiles excluding the interfacial regions of about 
4 a. As usual, the slip length was computed by linear extrapolation of the velocity profiles to 
the values = below the lower wall and Vx = U above the upper wall and then averaged. 



C. Shear viscosity and slip length 

In steady shear flow, the fluid viscosity is defined by the relation a^z = /^(t) 7, where 7 
denotes the shear rate and a^z is the shear stress through any plane parallel to the solid walls. 
In our simulations, the shear stress per unit area was computed at the liquid-solid interface 
by averaging the total force (in the shear flow direction) between the lower wall atoms and 
the fluid molecules. The dependence of polymer viscosity on shear rate is presented in Fig. [7] 
for selected values of the bending stiffness coefficient. For more flexible chains [kg ^ 2.0 e), 
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the gradual transition from the Newtonian to shear-thinning regimes is clearly observed 
in the accessible range of shear rates. The dashed line with the slope —0.37 is shown for 
reference in Fig.[71 indicating shear-thinning behavior of flexible chains {ko = 0.0 e and 



= 20) reported in previous studies 36|, |42| . The characteristic shear rate of the transition 



correlates well with the inverse relaxation time of polymer chains in the bulk region [see inset 
in Fig.|l](a)]. Not surprisingly, with increasing chain stiffness, the polymer viscosity at low 
shear rates increases, and the slope of the shear-thinning region becomes more steep due to 
the orientation of partially uncoiled chains in the shear flow. The apparent saturation of the 
viscosity at high shear rates is due to an increase in the fluid temperature near interfaces. 
The error bars are larger at low shear rates due to the enhanced statistical uncertainty in 
averaging velocity profiles and wall shear stress. 

FigureE] shows the dependence of slip length as a function of shear rate for the same flow 
conditions and values of the bending stiffness coefficient as in Fig. [71 All curves in Fig. [8] 
exhibit the same characteristic feature: a pronounced minimum at low shear rates and a 
steep increase at higher shear rates. In case of flexible chains, this behavior was analyzed 



previously [36|, l39|, |42| using Eq. ([T]). As illustrated in Fig.[Hl the slip length (for ko = 0.0 e) is 
nearly constant at low shear rates because of the extended Newtonian regime [in Fig. [7] and 
velocity-independent friction coefficient. With increasing shear rate, the relative competition 
between the shear-thinning viscosity and the dynamic friction coefficient in Eq. ([1]) leads to 
a minimum in the slip length, which is followed by a rapid increase at higher shear rates. 
Unexpectedly, increasing the chain stiffness produces larger slip lengths at low shear rates but 
smaller Lg at high shear rates. This trend can be understood by analyzing the effect of chain 
stiffness on the friction coefficient as a function of the slip velocity (see next subsection). 

As mentioned previously, the range of the upper wall speeds considered in the present 
study corresponds to anti-symmetric velocity profiles so that the slip velocity is the same at 
the lower and upper walls. It is expected that at higher upper wall speeds, the slip velocity 
at one of the solid walls will be much larger than the fluid thermal velocity producing slip 
lengths much larger than the channel height 38[. The investigation of the slip transition at 
very high shear rates is not the main focus of this paper; and, therefore, it was not studied 
in detail. 
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D. The dynamic friction coefficient 



In this subsection, we analyze the influence of chain stiffness on the friction coefficient 
at the liquid-solid interface as a function of the slip velocity. The results of previous MD 
studies have shown that the data for flexible polymer chains and weakly attractive crystalline 
surfaces can be well fitted by the following empirical equation: 

k/k* = + (12) 
where the parameter k* is the friction coefficient at small slip velocities and V* is the 



characteristic slip velocity of the transition to the nonlinear regime 



36 



39 



42|. It was 



demonstrated numerically that the friction coefficient k* is determined by the contact density 



and the in-plane structure factor of the first fiuid layer 



42l |. Furthermore, the characteristic 



slip velocity V* was found to correlate well with the diffusion rate of fiuid monomers over 



42| . The physical origin of 



the distance between nearest minima of the substrate potential 
the exponent —0.35 in Eq. f|T2|) is at present unclear. 

Although the friction coefficient can be readily computed from Eq. ([1]), the slight curvature 
in the velocity profiles near solid walls and the location of the liquid-solid interfaces used to 
compute the slip length [see Fig.|6]], introduce a small discrepancy between the definitions 
kiVs) = fi/Lg and k{Vi) = <Jxz/Vi, where Vg = Lg'j and Vi is the velocity of first fiuid layer. 
To eliminate this uncertainty, in the present study, the slip velocity was computed directly 
from the velocity profiles as follows: 

Vi= [ \^{z)p{z)dz I I \{z)dz, (13) 

J ZO J ZO 

where the limits of integration {zq = —11.54 a and zi = — 10.87 cr) define the width of the 
first peak in density profiles, which are shown for example in Fig. El 

The friction coefficient k{Vi) as a function of the slip velocity is plotted in Fig.|9] for 
several values of the bending stiffness coefficient. The important conclusion from the present 
results is that, with increasing chain stiffness, the friction coefficient at small slip velocities 
increases, and its decay rate at large slip velocities is independent of the chain stiffness. It 
can be further observed that for more flexible chains, kg ^ 1.0 e, the data are well described 
by the functional form given by Eq. ( IT2|) . However, as the chain stiffness increases, the data 
in Fig.|9] indicate qualitative changes in the velocity dependence of the friction coefficient. 



12 



i.e., the appearance of a pronounced maximum at small slip velocities. This non-monotonic 
behavior is related to the enhanced ordering of semifiexible chains near interfaces due to 
their orientation along the shear flow direction. This effect will be discussed in more detail 
in the next subsection. For the largest value of the stiffness coefficient, kg = 3.5 e, the error 
bars are relatively large at small slip velocities because the orientation of the extended chain 
segments in the first fluid layer is strongly influenced by the sixfold symmetry of the wall 
lattice and their relaxation dynamics is very slow [see Fig.|l](b)]. 

Nevertheless, some trends in the nonlinear rate dependence of the slip length presented in 
Fig.[H] can be understood from Eq. ([T]) and the data reported in Figs. [7] and IHl For example, 
the ratio of shear viscosity to the friction coefficient is smaller for stiffer chains with kg = 3.5 e 
at high shear rates, while the largest slip length at low shear rates is reported for chains 
with kg = 3.0 e. Also, the sharp decay of the slip length for the cases kg = 3.0 e and 3.5 e in 
Fig. [8] is related to the large negative slope of the polymer viscosity at low shear rates. As 
mentioned earlier, the nearly constant value of the slip length at low shear rates for fully 
flexible chains in Fig.|8]is due to the Newtonian viscosity and a wide linear regime of friction 
determined by the parameter V* in Eq. ( fT2i) . 



E. Fluid structure near solid walls 



The examples of the fluid density profiles shown in Fig. [5] demonstrate that the fluid 
density layering is most pronounced for the fluid monomers in contact with wall atoms. It is 
well known that, in addition to the fluid ordering perpendicular to the substrate, the periodic 
surface potential induces structure formation within the first fluid layer. The measure of the 
induced order is the in-plane static structure factor, which is defined as follows: 

1 1 2 

^ i=i 

where the sum is over A^^ fluid monomers in the layer and = [xj, yj) is the position vector of 
the j-th monomer. Depending on the strength of wall-fluid interactions and commensurabil- 
ity of liquid and solid structures, the structure factor typically contains a set of sharp peaks at 
the reciprocal lattice vectors, which are superimposed on several concentric rings characteris- 



tic of the liquid-like short range order {q]. In the past, several MD studies have demonstrated 
a strong correlation between the magnitude of the largest peak at the first reciprocal lattice 
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vector and the friction coefficient at liquid-solid interfaces |9|, m m m 128|, |35|, laei, 1391, |42|. 

We next plot the dependence of the normalized structure factor evaluated at the main 
reciprocal lattice vector Gi = (7.23 a^^, 0), the contact density, and the temperature of the 
ffist fluid layer in Fig-ITD] for three values of the bending stiffness coefficient. Similar to 
previous findings for flexible chains jsol, all three parameters in Fig. [TD] remain constant at 
small slip velocities, Vi < 0.1 ct/t, while the induced structure [S'(Gi)/S'(0) and pc] reduces 
and the fluid temperature increases at higher slip velocities. In sharp contrast, the structure 
factor for semiflexible chains exhibits a distinctive maximum at small slip velocities. Note 
that these changes in the structure factor are not reflected in the contact density, suggesting 
that they are mainly caused by the reorientation of chain segments in the first layer along 
the shear flow direction. 

In order to quantify this hypothesis, we examined the chain structure in contact with the 
substrate. Figure[TT] shows the average number of consecutive monomers per chain in the first 
fluid layer and their bond orientation with respect to the shear flow direction. Specifically, 
we computed the average value (cos^^), where 6 is the angle between the x axis and the 
three-dimensional bond vector connecting two consecutive monomers in the ffist fluid layer. 
In this definition, (cos^O) = 0.5 for the planar isotropic distribution, whereas (cos^O) = 1.0 
for the parallel arrangement of bond vectors along the x axis. The plots in Fig.lTTlreveal that, 
with increasing chain stiffness, the ffist fluid layer consists of more extended chain segments, 
which become preferentially aligned in the direction of shear flow. Notice that the orientation 
of flexible chain segments remain isotropic at small slip velocities, Vi < 0.1 a/r. Hence, the 
results in Figs.lTOl and ITT] indicate that the appearance of a maximum in 5'(Gi)/5'(0), which 
in turn affects the friction coefficient in Fig.[9l is due to the shear-induced alignment of 
semiflexible chain segments in the flrst fluid layer. 

We flnally summarize our data by plotting the inverse friction coefficient as a function 
of the combined variable S{0)/[S{Gi) pd in Fig.[12l It was previously shown for flexible 
polymer chains that in the linear regime [Vg < V* in Eq. (1121) ]. the friction coefficient can be 
described by a function of the variable S'(0)/[S'(Gi) pc] for a number of material parameters 
of the interface, such as fluid and wall densities, surface energy, chain length, and wall 
lattice type [42j]. The best flt to the MD data taken from 42| is indicated by the dashed 
line in Fig. [121 It can be observed that, for all values of the bending stiffness coefficient, 
the data points are distributed around the straight dashed line. In agreement with the 
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previous results [42!], at small slip velocities, the friction coefficient for more flexible chains, 
kg ^ 2.0 e, is well described by the master curve. The most noticeable difference between 
flexible and semiflexible chains in Fig.[T2]is the appearance of the hook-shaped curvature at 
small slip velocities, which is related to the local maximum in the structure factor discussed 
earlier. In other words, the same value of the product of structure factor and contact density, 
S'(0)/[S'(Gi) pc], corresponds to two different values of the friction coefficient, depending on 
the slip velocity and chain stiffness. 

In summary, the results in Fig.[T2]for semiflexible chains, kg ^ 2.0 e, confirm previous 
findings that the friction coefficient at small slip velocities is determined by the magnitude 
of the surface-induced peak in the structure factor and the contact density of the ffist fluid 



layer 



42| . The deviation from the master curve for more rigid chains, kg > 2.0 e, might be 



related to the slower relaxation dynamics of the chains in the interfacial region, similar to 
the trends found in dense polymer films at low shear rates 39 1. 



IV. CONCLUSIONS 



In this paper, we have presented results from extensive molecular dynamics simulations of 
thin polymer ffims confined by crystalline walls with weak surface energy. The computations 
were based on a coarse-grained bead-spring model of linear polymer chains with an additional 
bond angle potential that controls chain bending stiffness. The spatial configuration and 
local relaxation dynamics of polymer chains were characterized by the radius of gyration 
and the decay rate of the autocorrelation function of the ffist normal mode. We found that 
semiflexible chains near solid walls become more uncoiled and their relaxation dynamics is 
significantly slowed down. 

The most interesting result of the present study is the appearance of a distinct maximum 
in the velocity dependence of the friction coefficient due to the shear-induced alignment of 
semiflexible chain segments in the ffist fluid layer near solid walls. At small slip velocities, the 
orientation of more extended chain segments along the flow direction produces an enhanced 
ordering within the ffist fluid layer measured by the height of the main peak in the structure 
factor. This effect is absent for fully flexible chains since their segment orientation in the 
adjacent layer remains isotropic at small slip velocities. In addition, it was demonstrated 
that, with increasing slip velocity, the decay of the friction coefficient is independent of the 
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chain stiffness. 

Our simulation results indicate that the main features in the shear rate dependence of 
the slip length include a nearly constant value at low shear rates, a pronounced minimum 
at intermediate rates, and a rapid increase at high shear rates. These slip flow regimes are 
determined by the ratio of the rate-dependent polymer viscosity and the dynamic friction 
coefficient. Overall, we conclude that it is difficult to predict the net effect of chain stiff- 
ness on the slip length without performing numerical simulations; especially at low shear 
rates, where both polymer viscosity and friction coefficient increase with increasing bending 
rigidity. 
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FIG. 1: (Color online) A schematic view of the Couette flow configuration with slip at the lower 
and upper walls. Steady shear flow is generated by the upper wall moving with a constant velocity 
U in the x direction while the lower wall is at rest. The slip length and the slip velocity are related 
via Vs = 'yLs, where 7 is the shear rate computed from the slope of the velocity profile. 
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FIG. 2: (Color online) Instantaneous positions of fluid monomers (open blue circles) and fee wall 
atoms (filled gray circles) at equilibrium (i.e., both walls are at rest). Each monomer belongs to a 
polymer chain {N = 20) with the bending stiffness coefficient kg = 2.5 e. Seven chains are indicated 
by thick solid lines and filled black circles. The fluid monomer density is p = 0.91 cj~^ and the wall 
atom density is Pw — 1.40 fj 
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FIG. 3: (Color online) The averaged x, y, and z components of the radius of gyration Rgx (>), Rgy 
(v), Rgz (<l), and the total radius of gyration Rg (o) as a function of kg for polymer chains (a) in 
contact with the solid walls and (b) in the bulk region. The simulations were performed at the 
constant fluid density p = 0.91 while both walls were at rest. 
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FIG. 4: (Color online) The time autocorrelation function of the first normal mode Eq. (jlip for 
polymer chains (a) in the bulk region and (b) near the walls for several values of the bending 
stiffness coefficient. The inset shows the relaxation time of polymer chains in the bulk. 
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FIG. 5: (Color online) Averaged monomer density profiles near the lower stationary wall for the 
indicated values of the upper wall velocity U and p = 0.91 cj~'^. The bending stiffness coefficients 
are (a) ko = 0.0 e and (b) kg = 3.0 e. The left vertical axis at z = — 12.29 o" coincides with the fee 
lattice plane in contact with the polymer melt. The vertical dashed line at z = — 11.79 cr denotes 
the location of the liquid-solid interface. 



24 




FIG. 6: (Color online) Averaged normalized velocity profiles for the upper wall speeds (a) U = 
0.005 cr/r and (b) U = 0.5o"/r and bending stiffness coefficients kg = O.Oe (black lines), kg = 2.0e 
(red lines), and kg = 3.0 e (blue lines). The vertical axes coincide with the location of the fee lattice 
planes (at z/a = —12.29 and 9.73). The vertical dashed lines (at z/a = —11.79 and 9.23) indicate 
reference planes for computing the slip length. 
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FIG. 7: (Color online) Shear rate dependence of the polymer viscosity /x (in units of ercr^^) for 
selected values of the bending stiffness coefficient. The dashed line indicates a slope of —0.37. Solid 
curves are a guide for the eye. 
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FIG. 9: (Color online) Log-log plot of the friction coefficient k = axz/Vi (in units of ero"~^) as a 
function of the slip velocity of the first fluid layer Vi (in units of c/r) for the tabulated values of 
the bending stiffness coefficient. The dashed curve is the best fit to Eq. (112p with k* = 1.88 ero""'^ 
and V* = 0.2 (t/t. The solid curves are guides for the eye. 
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FIG. 10: (Color online) The normalized structure factor at the main reciprocal lattice vector 
Gi = (7.23 (T~^, 0) (a), contact density (b), and temperature (c) of the first fluid layer as a function 
of the slip velocity Vi (in units of ct/t). The values of the bending stiffness coefficient are kg = 0.0 e 
(o), ke = 2.0 e (<), and kg = 3.0 e (o). 
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FIG. 11: (Color online) The structure factor (a), number of consecutive monomers per chain (b), 
and bond orientation (c) in the first fluid layer as a function of the slip velocity. The values of 
the bending stiffness coefficient are kg = 0.0 e (o), kg = 2.0 e {<), and kg = 3.0 e (>). The data for 
S'(Gi)/S(0) are the same as in Fig.[TO](a). 
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FIG. 12: (Color online) Log-log plot of the inverse friction coefficient = Vijoxz (in units of 
cr^/er) as a function of Si^^j \S(Qi\) pc(y'^\ computed in the first fiuid layer. The values of the 
bending stiffness coefficient are tabulated in the inset. The dashed line y = 0.041 x^'^^ is taken 
from 421 and it is shown for reference. 
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